In this document, we will fit all models regarding this project. The detailed preregistered analysis plan can be found in our OSF’s project.
Let’s load libraries and data.
# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")
# Install/Load packages
pacman::p_load("rio",
"here",
"metafor",
"dplyr",
"bayesmeta",
"brms")
remotes::install_github("stan-dev/cmdstanr")
# Ensure that packages are installed according to versions in renv's lockfile
renv::restore()
# Import data
d = rio::import(here::here("data", "data.xlsx"))
Let’s calculate the crude odds ratio based on the extracted data.
# Calculate crude log odds ratios
d_logOR =
metafor::escalc(
measure = "OR", # log odds ratio,
# Tocilizumab
ai = trt_events,
n1i = trt_total,
# Control
ci = control_events,
n2i = control_total,
data = d
) %>%
as_tibble() %>%
# Calculate standard error
dplyr::mutate(sei = sqrt(vi),
# Set order to use "tocilizumab" as the Intercept in brms
treatment = factor(treatment,
levels = c("tocilizumab",
"sarilumab",
"baricitinib")))
# saveRDS(d_logOR,
# here("output", "data", "effect_sizes.Rds"))
First, we will fit a single meta-analysis model per drug (tocilizumab, sarilumab, or baricitinib vs. control).
All Bayesian meta-analysis random-effect models are the same, and is described as:
\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood} \\ \theta_i & \sim Normal(\mu, \tau^2)\\ \\ \mu & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors} \\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2) \\ \end{align*} \]
where \(y_i\) is the observed mean log odds ratio in study \(i\) for either tocilizumab, sarilumab, or baricitinib versus control treatment. We assume these effect sizes are normally distributed around the study-specific mean \(\theta_i\) along with a known sampling variance, represented by the observed \(\sigma_i^2\). We also assume \(\theta_i\) is drawn from a normal distribution where \(\mu\) is the average effect and \(\tau^2\) is the between-study heterogeneity.
## Formula
mf =
# https://bookdown.org/content/4857/horoscopes-insights.html#consider-using-the-0-intercept-syntax
formula(yi | se(sei) ~ 0 + Intercept + (1 | study))
## Priors
# Tau
informative = bayesmeta::TurnerEtAlPrior("all-cause mortality",
"pharmacological",
"placebo / control")
logmean = informative$parameters["tau", "mu"]
logsd = informative$parameters["tau", "sigma"]
# logmean
# -1.975
# logsd
# 0.67
priors_ma =
brms::prior(normal(0, 0.75), class = "b", coef = "Intercept") +
brms::prior(lognormal(-1.975, 0.67), class = "sd")
## Models
# Tocilizumab
ma_toci =
brms::brm(
data = d_logOR |> dplyr::filter(treatment == "tocilizumab"),
family = gaussian,
formula = mf,
prior = priors_ma,
sample_prior = TRUE,
control = list(adapt_delta = .95),
backend = "cmdstanr", # faster
cores = parallel::detectCores(),
chains = 4,
warmup = 2000,
iter = 4000,
seed = 123,
file = here::here("output", "fits", "ma_toci.Rds"),
file_refit = "on_change"
)
# Sarilumab
ma_sari =
brms::brm(
data = d_logOR |> dplyr::filter(treatment == "sarilumab"),
family = gaussian,
formula = mf,
prior = priors_ma,
sample_prior = TRUE,
control = list(adapt_delta = .95),
backend = "cmdstanr", # faster
cores = parallel::detectCores(),
chains = 4,
warmup = 2000,
iter = 4000,
seed = 123,
file = here::here("output", "fits", "ma_sari.Rds"),
file_refit = "on_change"
)
# Baricitinib
ma_bari =
brms::brm(
data = d_logOR |> dplyr::filter(treatment == "baricitinib", study != "RECOVERY Bari"),
family = gaussian,
formula = mf,
prior = priors_ma,
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
cores = parallel::detectCores(),
chains = 4,
warmup = 2000,
iter = 4000,
seed = 123,
file = here::here("output", "fits", "ma_bari.Rds"),
file_refit = "on_change"
)
ma_bari_sens =
brms::brm(
data = d_logOR |> dplyr::filter(treatment == "baricitinib", study != "RECOVERY Bari (No Toci)"),
family = gaussian,
formula = mf,
prior = priors_ma,
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
cores = parallel::detectCores(),
chains = 4,
warmup = 2000,
iter = 4000,
seed = 123,
file = here::here("output", "fits", "ma_bari_sens.Rds"),
file_refit = "on_change"
)
Now, we will fit Bayesian meta-regression models (one for tocilizumab vs. sarilumab data, another for tocilizumab vs. baricitinib).
The model can be described as:
\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood}\\ \theta_i & \sim Normal(\mu, \tau^2)\\ \mu &= \beta_0 + \beta_1 x_i\\ \\ \end{align*} \]
where \(y_i\) is the observed mean log odds ratio in study \(i\) for either tocilizumab, sarilumab, or baricitinib versus control treatment. We assume these effect sizes are normally distributed around the study-specific mean \(\theta_i\) along with a known sampling variance, represented by the observed \(\sigma_i^2\). We also assume \(\theta_i\) is drawn from a normal distribution where \(\mu\) is the average effect and \(\tau^2\) is the between-study heterogeneity. We will estimate \(\mu\) as a function of the estimated population effect size \(\beta_0\) and the moderator \(x\) multiplied by the coefficient \(\beta_1\). \(x\) is dummy-coded, where indicates that the observed \(y_i\) is an effect size estimate of compared to control, while indicates that \(y_i\) is an effect size estimate of the compared to control.
Here are the priors for this comparison. Of note, the “Optimistic for Sarilumab” and “Optimistic for Tocilizumab” models were the ones originally preregistered as main models. Yet they are presented as supplementary in the peer-reviewed article of this project.
\[ \begin{align*} \beta_0 & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors}\\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2)\\ \beta_{1[Skeptical]} & \sim \operatorname{Normal}(0, 0.354^2)\\ \beta_{1[Sarilumab]} & \sim \operatorname{Normal}(-0.049, 0.193^2)\\ \beta_{1[Tocilizumab]} & \sim \operatorname{Normal}(0.049, 0.193^2)\\ \beta_{1[Vague]} & \sim \operatorname{Normal}(0, 4^2) \end{align*} \]
## Means
# Figure S3 in https://doi.org/10.1101/2021.06.18.21259133;
# version posted June 22, 2021
reciprocal_remap_cap = 1/1.05
mean_log_sari = log(reciprocal_remap_cap)
mean_skeptical = log(1)
## SDs
# Assuming a mean = log(1), what is the SD that yields a distribution
# with 0.025 probability density below log(0.5)?
sari_sd_skeptical = (mean_skeptical - log(0.5))/qnorm(1-0.025)
# Assuming a mean = log(0.952), what is the SD that yields a distribution
# with 0.4 probability density above log(1)?
sari_sd_optimistic_sari = (log(1) - mean_log_sari)/qnorm(1-0.4)
# Assuming a mean = log(1/0.952) = -log(0.952), what is the SD that yields a
# distribution with 0.4 probability density below log(1)?
sari_sd_optimistic_toci = (-mean_log_sari - log(1))/qnorm(1-0.4)
sari_sd_vague = 4
# General model characteristics
## Formula
mf =
formula(yi | se(sei) ~ 0 + Intercept + treatment + (1 | study))
## Priors for beta_0 and tau
# tau:
# logmean
# -1.975
# logsd
# 0.67
priors =
brms::prior(normal(0, 0.75), class = "b", coef = "Intercept") +
brms::prior(lognormal(-1.975, 0.67), class = "sd", group = "study")
Fit models
## Skeptical model
# sari_sd_skeptical
# 0.353653
m_sari_skeptical =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentsarilumab"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_sari_skeptical.Rds"),
file_refit = "on_change"
)
## "Optimistic for Sarilumab" model
# mean_log_sari
# -0.04879016
# sari_sd_optimistic_sari
# 0.1925823
m_sari_optimistic_sari =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(-0.04879016, 0.1925823), class = "b", coef = "treatmentsarilumab"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_sari_optimistic_sari.Rds"),
file_refit = "on_change"
)
## "Optimistic for Tocilizumab" model
# -mean_log_sari
# 0.04879016
# sari_sd_optimistic_toci
# 0.1925823
m_sari_optimistic_toci =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0.04879016, 0.1925823), class = "b", coef = "treatmentsarilumab"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_sari_optimistic_toci.Rds"),
file_refit = "on_change"
)
## Vague model
m_sari_vague =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 4), class = "b", coef = "treatmentsarilumab"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_sari_vague.Rds"),
file_refit = "on_change"
)
Here are the priors for this comparison. Of note, the “Optimistic for Baricitinib” and “Optimistic for Tocilizumab” models were the ones originally preregistered as main models. Yet they are presented as supplementary in the peer-reviewed article of this project.
\[ \begin{align*} \beta_0 & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors}\\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2)\\ \beta_{1[Skeptical]} & \sim \operatorname{Normal}(0, 0.354^2)\\ \beta_{1[Baricitinib]} & \sim \operatorname{Normal}(-0.105, 0.416^2)\\ \beta_{1[Tocilizumab]} & \sim \operatorname{Normal}(0.105, 0.416^2)\\ \beta_{1[Vague]} & \sim \operatorname{Normal}(0, 4^2) \end{align*} \]
## Means
mean_log_bari = log(0.9)
mean_skeptical = log(1)
## SDs
# Assuming a mean = log(1), what is the SD that yields a distribution
# with 0.025 probability density below log(0.5)?
bari_sd_skeptical = (mean_skeptical - log(0.5))/qnorm(1-0.025)
# Assuming a mean = log(0.9), what is the SD that yields a distribution
# with 0.4 probability density above log(1)?
bari_sd_optimistic_bari = (log(1) - mean_log_bari)/qnorm(1-0.4)
# Assuming a mean = log(1/0.9) = -log(0.9), what is the SD that yields a
# distribution with 0.4 probability density below log(1)?
bari_sd_optimistic_toci = (-mean_log_bari - log(1))/qnorm(1-0.4)
bari_sd_vague = 4
## Skeptical model
# bari_sd_skeptical
# 0.353653
m_bari_skeptical =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_skeptical.Rds"),
file_refit = "on_change"
)
## "Optimistic for Baricitinib" model
# mean_log_bari
# -0.1053605
# bari_sd_optimistic_bari
# 0.4158742
m_bari_optimistic_bari =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(-0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_optimistic_bari.Rds"),
file_refit = "on_change"
)
## "Optimistic for Tocilizumab" model
# -mean_log_bari
# 0.1053605
# bari_sd_optimistic_toci
# 0.4158742
m_bari_optimistic_toci =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_optimistic_toci.Rds"),
file_refit = "on_change"
)
## Vague model
m_bari_vague =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 4), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_vague.Rds"),
file_refit = "on_change"
)
Here we will use the post-hoc priors based on Karampitsakos et al.’s results for the “Optimistic for Baricitinib” and “Optimistic for Tocilizumab” models.
# Data provided by Karampitsakos et al. through email
# Calculate crude log odds ratio
bari_toci_data =
metafor::escalc(
measure = "OR", # log odds ratio
# Baricitinib
ai = 40,
n1i = 125,
# Tocilizumab
ci = 50,
n2i = 126
) %>%
as_tibble() %>%
# Calculate standard error
dplyr::mutate(sei = sqrt(vi))
## Mean
mean_log_bari = bari_toci_data$yi[1]
## SDs
bari_sd_optimistic_bari = bari_toci_data$sei[1]
bari_sd_optimistic_toci = bari_sd_optimistic_bari
## Skeptical model
# bari_sd_skeptical
# 0.353653
m_bari_skeptical_sens =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_skeptical_sens.Rds"),
file_refit = "on_change"
)
## "Optimistic for Baricitinib" model
# mean_log_bari
# -0.3350615
# bari_sd_optimistic_bari
# 0.2644288
m_bari_optimistic_bari_direct_sens =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(-0.3350615, 0.2644288), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_optimistic_bari_direct_sens.Rds"),
file_refit = "on_change"
)
## "Optimistic for Tocilizumab" model
# -mean_log_bari
# 0.3350615
# bari_sd_optimistic_toci
# 0.2644288
m_bari_optimistic_toci_direct_sens =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0.3350615, 0.2644288), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_optimistic_toci_direct_sens.Rds"),
file_refit = "on_change"
)
## Vague model
m_bari_vague_sens =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 4), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_vague_sens.Rds"),
file_refit = "on_change"
)
# Load function
source(here::here("functions", "diag_plot.R"))
## * The library is already synchronized with the lockfile.
Tocilizumab model
diag_plot(model = ma_toci,
pars_list = c("b_Intercept", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(ma_toci,
type = "dens_overlay",
ndraws = 50)
brms::pp_check(ma_toci,
type = "ecdf_overlay",
ndraws = 50)
Sarilumab model
diag_plot(model = ma_sari,
pars_list = c("b_Intercept", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(ma_sari,
type = "dens_overlay",
ndraws = 50)
brms::pp_check(ma_sari,
type = "ecdf_overlay",
ndraws = 50)
Baricitinib (sensitivity) models
diag_plot(model = ma_bari,
pars_list = c("b_Intercept", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(ma_bari,
type = "dens_overlay",
ndraws = 50) + coord_cartesian(y = c(0, 5))
brms::pp_check(ma_bari,
type = "ecdf_overlay",
ndraws = 50)
Baricitinib models
diag_plot(model = ma_bari_sens,
pars_list = c("b_Intercept", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(ma_bari_sens,
type = "dens_overlay",
ndraws = 50) + coord_cartesian(y = c(0, 5))
brms::pp_check(ma_bari_sens,
type = "ecdf_overlay",
ndraws = 50)
“Skeptical” model
diag_plot(model = m_sari_skeptical,
pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_sari_skeptical,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50)
brms::pp_check(m_sari_skeptical,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Optimistic for Sarilumab” model
diag_plot(model = m_sari_optimistic_sari,
pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_sari_optimistic_sari,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50)
brms::pp_check(m_sari_optimistic_sari,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Optimistic for Tocilizumab” model
diag_plot(model = m_sari_optimistic_toci,
pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_sari_optimistic_toci,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50)
brms::pp_check(m_sari_optimistic_toci,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Vague” model
diag_plot(model = m_sari_vague,
pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_sari_vague,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50)
brms::pp_check(m_sari_vague,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Skeptical” models
diag_plot(model = m_bari_skeptical,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_skeptical,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_skeptical,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
diag_plot(model = m_bari_skeptical_sens,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_skeptical_sens,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_skeptical_sens,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Optimistic for Baricitinib” models
diag_plot(model = m_bari_optimistic_bari,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_optimistic_bari,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_optimistic_bari,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
diag_plot(model = m_bari_optimistic_bari_direct_sens,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_optimistic_bari_direct_sens,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_optimistic_bari_direct_sens,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Optimistic for Tocilizumab” models
diag_plot(model = m_bari_optimistic_toci,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_optimistic_toci,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_optimistic_toci,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Optimistic for Tocilizumab” model
diag_plot(model = m_bari_optimistic_toci_direct_sens,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_optimistic_toci_direct_sens,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_optimistic_toci_direct_sens,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Vague” models
diag_plot(model = m_bari_vague,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_vague,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_vague,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
diag_plot(model = m_bari_vague_sens,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_vague_sens,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_vague_sens,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)